# carregando pacotes
library(tidyverse)
library(factoextra)
library(FactoMineR)
library(lme4)
library(stargazer)
library(cowplot)
library(ggpubr)
library(broom)
library(rgdal)
library(mapproj)
library(stringi)
library(sjPlot)
library(sjmisc)
library(psych)
library(scales)
library(ltm)
library(mirt)
library(ggcorrplot)
library(GGally)
library(texreg)
library(radiant.model)
library(dotwhisker)
library(Zelig)
library(arm)
library(car)
library(reghelper)

setwd('C:/Users/borba/OneDrive/�rea de Trabalho/2019/dissertation')

##### DADOS #####

# ess
ess <- read.csv('./Analysis data/ess2.csv')

attach(ess)

# clea - somente partidos de direita radical
clea <- read.csv('./Analysis data/clea_rrp.csv')

# bundeswahlleiter2017
bw <- read.csv('./Original data/btw17.csv', sep = ';')

# georgiadou et at
farright_data<-read.dta(file="GRR_dataMapsR.dta")

##### criando fun��es #####
# gr�fico de densidade
dplot <- function(VB, VW, Numero, Titulo) {
        ggplot() +
                geom_density(data = essgen, aes(x = VB, fill = "black"), alpha = .3,
                             adjust = 3.5) +
                geom_density(data = essrrp, aes(x = VW, fill = "grey"), alpha = .3, adjust = 1.5) +
                scale_x_continuous(name = "", breaks = 0 : 10, limits = c(0,10), expand = c(.01,.01)) +
                scale_y_continuous(name = "", breaks =  NULL, labels = NULL, limits = c(0, .5)) +
                scale_fill_identity(name = "", labels = c("Parties", "Voters")) +
                ggtitle(Titulo) +
                theme_classic(base_size = 12) +
                theme(legend.position = c(.3, .7), legend.background = element_blank(), 
                      legend.text = element_text(size = 7), 
                      plot.title = element_text(size = 10, hjust = 0.5, face = "bold"),
                      plot.caption = element_text(hjust = .5),
                      axis.text.x = element_text(size = 7))
}

# gr�fico de densidade com eixo y menor
dplotb <- function(VB, VW, Numero, Titulo) {
        ggplot() +
                geom_density(data = essgen, aes(x = VB, fill = "black"), alpha = .3,
                             adjust = 7) +
                geom_density(data = essrrp, aes(x = VW, fill = "grey"), alpha = .3, adjust = 5) +
                scale_x_continuous(name = "", breaks = 0 : 4, limits = c(0,4), expand = c(.01,.01)) +
                scale_y_continuous(name = "", breaks =  NULL, labels = NULL, limits = c(0, .7)) +
                scale_fill_identity(name = "", labels = c("Parties", "Voters")) +
                ggtitle(Titulo) +
                theme_classic(base_size = 12) +
                theme(legend.position = c(.3, .7), legend.background = element_blank(), 
                      legend.text = element_text(size = 7), 
                      plot.title = element_text(size = 10, hjust = 0.5, face = "bold"),
                      plot.caption = element_text(hjust = .5),
                      axis.text.x = element_text(size = 7))
}

# multiplot 
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
        library(grid)
        
        # Make a list from the ... arguments and plotlist
        plots <- c(list(...), plotlist)
        
        numPlots = length(plots)
        
        # If layout is NULL, then use 'cols' to determine layout
        if (is.null(layout)) {
                # Make the panel
                # ncol: Number of columns of plots
                # nrow: Number of rows needed, calculated from # of cols
                layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                                 ncol = cols, nrow = ceiling(numPlots/cols))
        }
        
        if (numPlots==1) {
                print(plots[[1]])
                
        } else {
                # Set up the page
                grid.newpage()
                pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
                
                # Make each plot, in the correct location
                for (i in 1:numPlots) {
                        # Get the i,j matrix positions of the regions that contain this subplot
                        matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
                        
                        print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                                        layout.pos.col = matchidx$col))
                }
        }
}

##### HISTOGRAMAS E GR�FICOS DE BARRA #####

props <- ddply(.data = ess, .var = "cntry", .fun = function(x) {
        data.frame(n = nrow(x),
                   imwbcnt = nrow(subset(x, jc13 == 'Yes, a military take-over of the state would be justified')),
                   imwbcnt.prop = nrow(subset(x, jc13 == 'Yes, a military take-over of the state would be justified')) / nrow(x)
        )
}
)

prop.table(ess$imwbcnt)
ddpl

props <- ess %>% group_by(cntry) %>% data.frame(n = nrow(ess),
                                                imwbcnt.prop = nrow(ess$imwbcnt)/nrow(ess))

# gr�fico de barra da média de votos por partido
hist1 <- ggplot(data = clea, aes(x = reorder(pty_n, +pvs1), y = pvs1)) +
        stat_summary(fun.y = mean, geom = "bar", aes(group = 1)) +
        coord_flip() +
        theme_classic() +
        labs(x = "Partido", y = "Média de votos (%)") +
        theme(plot.title = element_text(size = 12, hjust = .5, face = "bold"), 
              plot.caption = element_text(hjust = .5)) +
        geom_hline(aes(yintercept = mean(clea$pvs1)), linetype = 2)

ggsave("./Command files/plots/hist_partymean.png", plot = hist1, device = "png", width = 5, height = 5)

# histograma da vari�vel dependente
hist2 <- ggplot(data = ess, aes(x = as.factor(vote_rrp))) +
        geom_bar(width = .25) +
        theme_sjplot2() +
        labs(x = '', y = '') +
        scale_y_continuous(breaks = seq(0, 124000, 10000))

ggsave('./Command files/plots/hist_vote_rrp.png', plot = hist2, device = 'png',
       width = 4, height = 3)

# histograma da vari�vel dependente por pa�s
describeBy(vote_rrp, cntry)

hist3 <- ggplot(data = ess, aes(x = vote_rrp, fill = cntry)) +
        geom_bar(width = .9, position = 'dodge2') +
        theme_sjplot2() +
        labs(x = '', y = '') +
        scale_x_continuous(breaks = c(0, 1)) +
        scale_y_continuous(breaks = seq(0, 25000, 2500)) +
        scale_fill_discrete(name = 'Pa�s')

ggsave('./Command files/plots/hist_vote_rrp_ctry.png', plot = hist3, device = 'png',
       width = 5, height = 4)

# invertendo valores da variável
ess$imwbcnt2 <- max(imwbcnt, na.rm = T) - imwbcnt

# histograma da variável independente
hist4 <- ggplot(data = ess, aes(x = imwbcnt2)) +
        geom_histogram(position = 'dodge2', binwidth = 1, colour = 'white') +
        theme_sjplot2() +
        labs(x = '', y = '') +
        scale_x_continuous(breaks = seq(0, 10, 1))

ggsave('./Command files/plots/hist_imwbcnt.png', plot = hist4, device = 'png',
       width = 5, height = 4)


# gráfico de barra da variável independente por pa�s
hist5 <- ggplot(data = ess, aes(x = imwbcnt2, fill = cntry)) +
        geom_bar(position = 'dodge2') +
        theme_sjplot2() +
        labs(x = '', y = '') +
        scale_fill_discrete(name = 'Pa�s') +
        scale_x_continuous(breaks = seq(0, 10, 1)) +
        scale_y_continuous(breaks = seq(0, 8000, 1000))

ggsave('./Command files/plots/hist_imwbcnt_cntry.png', plot = hist5, device = 'png',
       width = 5, height = 4)

# histograma das vari�veis para a cria��o do �ndice
hist6 <- ggplot(data = ess, aes(x = imsmetn)) +
        geom_histogram(position = 'dodge2', binwidth = 1, colour = 'white') +
        theme_sjplot2() +
        labs(x = '', y = '') +
        scale_x_continuous(breaks = seq(1, 10, 1)) +
        scale_y_continuous(breaks = seq(0, 80000, 10000))

# invertendo valores
ess$imbgeco2 <- max(imbgeco, na.rm = T) - imbgeco
ess$imueclt2 <- max(imueclt, na.rm = T) - imueclt

# fun��oo para os histogramas
h <- function(x){
        ggplot(data = ess, aes(x = x)) +
                geom_histogram(position = 'dodge2', binwidth = 1, colour = 'white') +
                theme_sjplot2() +
                labs(x = '', y = '') +
                scale_x_continuous(breaks = seq(0, 10, 1))
}

hist7 <- h(imdfetn) 
hist8 <- h(impcntr)
hist9 <- h(ess$imbgeco2)
hist10 <- h(ess$imueclt2)

save <- function(name, plot){
        ggsave(filename = name, plot = plot, device = 'png', width = 5, height = 4)
}

save('./Command files/plots/imsmetn.png', hist6)
save('./Command files/plots/imdfetn.png', hist7)
save('./Command files/plots/impcntr.png', hist8)
save('./Command files/plots/imbgeco.png', hist9)
save('./Command files/plots/imuelct.png', hist10)

# histograma das variáveis de controle
# variáveis de satisfação
hist11 <- h(stflife)
hist12 <- h(stfdem)
hist13 <- h(stfeco)
hist14 <- h(stfgov)

# salvando histogramas
save('./Command files/plots/stflife.png', hist11)
save('./Command files/plots/stfdem.png', hist12)
save('./Command files/plots/stfeco.png', hist13)
save('./Command files/plots/stfgov.png', hist14)

# número de imigrantes no país
hist15 <- h(noimbro) +
        scale_x_continuous(breaks = seq(0, 100, 10))

save('./Command files/plots/noimbro.png', hist15)

ggplot(data = ess) +
        geom_density(aes(x = noimbro, fill = 'red'), alpha = .3) +
        geom_density(aes(x = noimbro2, fill = 'gray'), alpha = .3)
        

# proximidade com partidos
hist16 <- ggplot(data = ess, aes(x = partisan)) +
        geom_histogram(position = 'dodge2', binwidth = .5) +
        scale_x_continuous(breaks = c(0, 1)) +
        theme_sjplot2() +
        labs(x = '', y = '')

save('./Command files/plots/partisan')

summary(agea)

##### BOXPLOTS #####
# boxplot da distribuição dos votos dos partidos da direita radical por país
box1 <- ggplot(data = clea, aes(x = reorder(ctr_n, +pvs1), y = pvs1)) +
        geom_boxplot() +
        stat_summary(fun.y = mean, geom = "point", colour = "red", size = 1.5) +
        coord_flip() +
        theme_classic() +
        labs(x = "País", y = "Porcentagem de votos") +
        theme(plot.title = element_text(size = 12, hjust = .5, face = "bold"),
              plot.caption = element_text(hjust = .5)) +
        scale_y_continuous(breaks = seq(.0, .7, .05))

ggsave("./Command files/plots/box_rrpp_by_ctr.png", plot = box1, device = "png", width = 6, height = 4)

# criando variável dummy que sinaliza um land do leste alemão
bw$east <- as.factor(ifelse(bw$Land %in% c(13, 12, 15, 14, 16), 1, 0))

bw$AfD_prop <- bw$AfD / bw$Gültige

# boxplot do desempenho da AfD nas eleições de 2017
box2 <- ggplot(bw, aes(x = factor(Land), y = AfD, fill = east)) +
        geom_boxplot() +
        theme_minimal() +
        theme(plot.title = element_text(size = 12, face = "bold", hjust = .5),
              plot.caption = element_text(hjust = .5)) +
        labs(x = "Região", y = "Número de votos AfD")

# plotando boxplot
box3 <- ggplot(bw, aes(x = factor(Land), y = AfD_prop * 100, fill = east)) +
        geom_boxplot() +
        theme_minimal() +
        theme(plot.title = element_text(size = 12, face = "bold", hjust = .5),
              plot.caption = element_text(hjust = .5)) +
        labs(x = "Região", y = "Votos AfD (%)") +
        scale_fill_discrete(name = 'Leste') +
        scale_y_continuous(breaks = seq(0, 60, 10))

ggsave('./plots/boxplot_afd17.png', plot = box3, device = 'png', width = 6, height = 3.5)

# boxplot com os dados de georgiadou et al
farright_data$country <- stri_sub(farright_data$id, from = 0, to = 2)

box4 <- ggplot(data = farright_data, aes(x = reorder(country, -prr_average), y = prr_average)) +
        geom_boxplot() +
        stat_summary(fun.y = mean, geom = "point", colour = "red", size = 1.5) +
        coord_flip() +
        theme_classic() +
        labs(x = "País", y = "Porcentagem de votos") +
        theme(plot.title = element_text(size = 12, hjust = .5, face = "bold"),
              plot.caption = element_text(hjust = .5))

# gráfico de coeficientes - teste-t entre leste e oeste alemão
# criando função para cálculo de intervalo de confiança
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=TRUE,
                      conf.interval=.95, .drop=TRUE) {
        library(plyr)
        
        # New version of length which can handle NA's: if na.rm==T, don't count them
        length2 <- function (x, na.rm=FALSE) {
                if (na.rm) sum(!is.na(x))
                else       length(x)
        }
        
        # This does the summary. For each group's data frame, return a vector with
        # N, mean, and sd
        datac <- ddply(data, groupvars, .drop=.drop,
                       .fun = function(xx, col) {
                               c(N    = length2(xx[[col]], na.rm=na.rm),
                                 mean = mean   (xx[[col]], na.rm=na.rm),
                                 sd   = sd     (xx[[col]], na.rm=na.rm)
                               )
                       },
                       measurevar
        )
        
        # Rename the "mean" column    
        datac <- rename(datac, c("mean" = measurevar))
        
        datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
        
        # Confidence interval multiplier for standard error
        # Calculate t-statistic for confidence interval: 
        # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
        ciMult <- qt(conf.interval/2 + .5, datac$N-1)
        datac$ci <- datac$se * ciMult
        
        return(datac)
}



# plotando gráfico
ggplot(data = t_test_data, aes(y = estimate, x = as.factor(east))) +
        geom_point() +
        geom_errorbar(aes(ymin = estimate - ci, ymax = estimate + ci), width = .1) +
        labs(y = "Média de votos AfD (%)", x = "Leste") +
        theme_classic() +
        theme(plot.title = element_text(face = "bold", hjust = .5, size = 12),
              plot.caption = element_text(hjust = .5)) +
        scale_y_continuous(breaks = seq(0, 30, 5))




##### GRÁFICOS DE LINHA #####
# calculando média de votos dos partidos por ano
plot <- clea %>% group_by(pty_n, yr, ctr_n) %>% summarise_at(vars(pvs1), 
                                                funs(mean(., na.rm = T))) %>%
        arrange(desc(pvs1))

# visualizando os cinco partidos com maior média de votos
agg <- aggregate(plot, by = list(plot$pty_n), FUN = mean, na.rm = T) %>% 
        arrange(desc(pvs1))
head(agg, n = 5)

# agregando média por país
agg_ctry <- aggregate(clea, by = list(clea$ctr_n), FUN = mean, na.rm = T)

# plotando a média de votos dos partidos por ano
line2 <- ggplot(plot %>% filter(pty_n %in% c('AN', 'FPÖ', 'FN', 'PVV', 'LN'))) +
         geom_line(aes(x = yr, y = pvs1*100, colour = pty_n)) +
         theme_minimal() +
         labs(x = '', y = 'Votos (%)') +
         scale_color_discrete(name = 'Partidos') +
         scale_x_continuous(limits = c(1955, 2018), breaks = seq(1955, 2018, 10), 
                            labels = seq(1955, 2018, 10))

ggsave("timeseries2_rrpp.png", plot = line2, device = "png", width = 6, height = 3.5)

# gráfico de linha do desempenho dos maiores cinco partidos
line1 <- ggplot() +
        stat_summary(data = clea_rrp[clea_rrp$ctr_n == "Finland",], aes(x = yr, y = pvs1),
                     fun.y = mean, geom = "line", color = "blue") +
        stat_summary(data = clea_rrp[clea_rrp$ctr_n == "Italy" & clea_rrp$pty_n %in% c("LN", "MSI"),], aes(x = yr, y = pvs1),
                     fun.y = mean, geom = "line", color = "green") +
        stat_summary(data = clea_rrp[clea_rrp$ctr_n == "Austria" & clea_rrp$pty_n == "FPÖ",],
                     aes(x = yr, y = pvs1), fun.y = mean, geom = "line", color = "red") +
        stat_summary(data = clea_rrp[clea_rrp$ctr_n == "France" & clea_rrp$pty_n == "FN",],
                     aes(x = yr, y = pvs1), fun.y = mean, geom = "line", color = "lightblue") +
        stat_summary(data = clea_rrp[clea_rrp$ctr_n == "Netherlands",],
                     aes(x = yr, y = pvs1), fun.y = mean, geom = "line", color = "orange") +
        theme_classic() +
        labs(x = "Ano da eleição", y = "Média de votos (%)", caption = "Figura 3
             Elaboração do autor com dados do CLEA") +
        ggtitle("Média de votos na direita radical (cinco primeiros países)") +
        theme(plot.title = element_text(size = 12, hjust = .5, face = "bold"),
              plot.caption = element_text(hjust = .5)) +
        scale_x_continuous(limits = c(1950, 2017), breaks = seq(1950, 2017, 5)) +
        geom_vline(xintercept = 2007, linetype = 2, color = "grey 60")

ggsave("timeseries_rrpp.png", plot = line5, device = "png", width = 6, height = 3.5)

# série histórica - média geral
line3 <- ggplot(data = clea, aes(x = yr, y = pvs1)) +
        stat_summary(fun.y = mean, geom = 'line')

##### MAPA #####

# shape
shapefile_w <- readOGR("../Original data/shapefiles/Europe.shp")

shapefile_w@data

# convertendo o shapefile para dataframe 
shapefile_df <- fortify(shapefile_w)

shapefile_data <- fortify(shapefile_w@data)

# criando variável id para fazer o merge com os países
shapefile_data$id <- row.names(shapefile_data) 

# merge com os países
shapefile_df <- full_join(shapefile_df, shapefile_data, by="id")

# criando variável para fazer o merge com os dados do clea
agg_ctry$NAME <- agg_ctry$Group.1 %>% as.character()
agg_ctry$NAME[agg_ctry$NAME == 'UK'] <- 'United Kingdom'

# transformando variável país para caracter
shapefile_df$NAME <- shapefile_df$NAME %>% as.character()

# selecionando países
print(agg_ctry$NAME) %>% as.character()
shapefile_df <- shapefile_df %>% 
        filter(!NAME %in% c('Faeroe Islands (Denmark)', 'Gibraltar (UK)',
                           'Guernsey (UK)', 'Iceland', 'Isle of Man (UK)',
                           "Jan Mayen (Norway)", 'Jersey (UK)', 'Svalbard (Norway)',
                           'Russia', 'Turkey', 'Georgia', 'Armenia', 
                           'Azerbaijan'))

# mergindo dados
shapefile_df <- left_join(shapefile_df, agg_ctry %>% select(NAME, pvs1), 
                          by = 'NAME')

# plotando mapa
map1 <-  ggplot() + geom_polygon(data = shapefile_df,
                                aes(x = long, y = lat, group = group, 
                                    fill = pvs1 * 100),
                                colour = "grey", size = .1) +
        theme_void(base_size = 11) +
        scale_fill_viridis(option = 'magma',
                           name = 'Votos (%)', direction = -1,
                           guide = guide_colorbar(
                                   direction = "horizontal",
                                   barheight = unit(2, units = "mm"),
                                   barwidth = unit(50, units = "mm"),
                                   draw.ulim = F,
                                   title.position = 'top',
                                   title.hjust = 0.5,
                                   label.hjust = 0.5
                           )) +
        coord_map() +
        theme(legend.position = "bottom",
              plot.title = element_text(size = 12, face = 'bold', hjust = .5),
              plot.caption = element_text(hjust = .5))

ggsave('../Command files/plots/map_rr_europe.png', plot = map1, device = 'png', width = 6, height = 5)

# mapa com os dados de georgiadou et al
EU_NUTS <- readOGR(dsn = "../replication_map", layer = "NUTS_RG_60M_2010")

### Change projection to Google-like Map
proj4string(EU_NUTS)

proj4string(EU_NUTS) <- CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs")

EU_NUTS <- spTransform(EU_NUTS, CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"))

plot(EU_NUTS,xlim = c(300000,1900100), ylim=c(3500000,11200000))
##EU_NUTS <- spTransform(x, CRS("+init=epsg:4238"))


farright_data$NUTS_ID<-farright_data$nuts_2   #create common ID variable for merging spatial and political datasets

farright_data<-farright_data[which(farright_data$year=="2013"),] #get rid of multiple years (we are plotting averages by nuts 2 through time)

EU_NUTS2013_fortified <- fortify(EU_NUTS, region = "NUTS_ID") ##create fortified databse for plotting in ggplot

farright_data$id<-as.character(farright_data$nuts_2)  ## create character id in political variable for merging with geographical database

# criando média dos votos dos dois grupos de partidos
farright_data$rad_right <- (farright_data$prr_average + farright_data$er_average) / 2

# mergindo dados geográficos com dados eleitorais
map_data <-left_join(EU_NUTS2013_fortified,farright_data, by = c("id", "id"))  ## merge political and geographical databases

map_data$nuts2 <- as.character(map_data$nuts_2)

map_data2 <- map_data[nchar(map_data$nuts2)==4,]
map_data3 <- map_data[nchar(map_data$nuts2)==4,]

# plotando mapa
q <- ggplot() +
        geom_polygon(data = map_data2 %>% filter(!id %in% c('FR90', 'FR91', 'FR92',
                                                            'FR93', 'FR94')),
                     aes(fill = prr_average * 100, x = long, y = lat, group = group)) +
        theme_void() +
        coord_map() +
        scale_fill_viridis(option = 'magma', na.value = 'grey',
                           name = 'Votos (%)', direction = -1,
                           guide = guide_colorbar(
                                   direction = "horizontal",
                                   barheight = unit(2, units = "mm"),
                                   barwidth = unit(50, units = "mm"),
                                   draw.ulim = F,
                                   title.position = 'top',
                                   title.hjust = 0.5,
                                   label.hjust = 0.5
                           )) +
        theme(legend.position = 'bottom')

# exportando
ggsave('../Command files/plots/maprr.png', plot = q, device = 'png', width = 6, height = 5)

##### gráficos de densidade #####
# gráficos de densidade para o posicionamento do eleitorado e dos partidos
# invertendo valores das variáveis de percepção sobre imigração (valores altos
# indicam opinião contrária aos imigrantes)
ess$imwbcnt2 <- max(ess$imwbcnt, na.rm = T) - ess$imwbcnt
ess$imueclt2 <- max(ess$imueclt, na.rm = T) - ess$imueclt
ess$imbgeco2 <- max(ess$imbgeco, na.rm = T) - ess$imbgeco

table(ess$imbgeco)

# criando banco de dados somente com observações de eleitores de RRPPs
essrrp <- ess[ess$vote_rrp == 1,]

# dados para eleitorado geral
essgen <- ess[ess$vote_rrp == 0,]

# criando banco de dados para o eleitorado geral (exceto eleitores de RRPPs)
# density plot 1 - posicionamento dos eleitores e partidos num eixo esquerda-direita
dplot1 <- dplot(essgen$lrscale, essrrp$lrscale, 1, 
                "Auto posicionamento escala esquerda-direita") +
        scale_fill_identity(name = "", labels = c("Eleitorado", "Eleitor RRPP"), 
                            guide = "legend")

# density plot 2 - imigrantes fazem do país um lugar melhor ou pior para se viver
dplot2 <- dplot(essgen$imwbcnt, essrrp$imwbcnt, 2, "Imigrantes fazem do país um lugar pior")

# density plot 3 - sentimento do eleitorado sobre o impacto cultural do imigrante
dplot3 <- dplot(essgen$imueclt, essrrp$imueclt, 3, "Imigrantes prejudicam a cultura do país")

# density plot 4 - importância de seguir costumes e tradições
dplot4 <- dplotb(essgen$imdfetn, essrrp$imdfetn, 4, "Restrição de imigrantes de grupos étnicos diferentes")

# exportando gráficos
dplots <- ggarrange(dplot1, dplot2, dplot3, dplot4, 
                    ncol = 2, nrow = 2)
ggsave('./plots/dplots.png', plot = dplots, device = 'png', 
       width = 8, height = 6)

##### TABELAS #####

##### Testes-t
# teste-t para as variáveis 
tt_lrscale <- t.test(essrrp$lrscale, essgen$lrscale) 
tt_imw <- t.test(essgen$imwbcnt, essrrp$imwbcnt)
tt_imu <- t.test(essgen$imueclt, essrrp$imueclt)
tt_imd <- t.test(essrrp$imdfetn, essgen$imdfetn)

# convertendo os objetos para data dataframe
tt_lrscale <- tt_lrscale %>% tidy()
tt_imw <- tt_imw %>% tidy()
tt_imu <- tt_imu %>% tidy()
tt_imd <- tt_imd %>% tidy()

# juntando os testes t
t <- bind_rows(tt_lrscale, tt_imw, tt_imu, tt_imd)

# teste-t da média de votos da AfD no leste e oeste alemão
t1 <- t.test(bw$AfD_prop ~ bw$east)

t1_tidy <- tidy(t1)

# construindo bancos de dados com os resultados do teste-t
t_test_data <- data.frame(estimate = c(11.24, 23.53), ci = rep(.1241), east = c(0,1))

# função para renomear linhas e colunas da tabela
stargazer_htest = function (data, ...) {
        summary = data.frame(`T-statistic` = data$statistic,
                             DF = data$parameter,
                             `P value` = data$p.value,
                             `Alternative hypothesis` = data$alternative,
                             check.names = FALSE)
        stargazer(summary, flip = TRUE, summary = FALSE,
                  notes = paste(data$method, data$data.name, sep = ': '), ...)
}

# stargazer t-test
stargazer_htest(t, column.labels = c("lrscale", "imwbcnt", "imueclt", "imdfetn"), 
                out = "t-test.tex")

stargazer_htest(t1_tidy, out = 't-test_afd.tex')

##### MODELOS IRT #####
#ess$imsmetn <- ess$imsmetn -1
#ess$imdfetn <- ess$imdfetn -1
#ess$impcntr <- ess$impcntr -1
#ess$imbgeco <- ess$imbgeco +1
#ess$imueclt <- ess$imueclt +1
#ess$imwbcnt <- ess$imwbcnt +1

# removendo NA's das variáveis selecionadas
#ess <- ess %>% drop_na(imsmetn, imdfetn, impcntr, imbgeco2, imueclt2, imwbcnt2)

# estimando graded response model para as variáveis sobre imigração
irt <- grm(ess[, c("imsmetn", "imdfetn", "impcntr", "imbgeco2", "imueclt2", 
                   "imwbcnt2")])

cronbach.alpha(ess[, c("imsmetn", "imdfetn", "impcntr", "imbgeco2", "imueclt2", 
                       "imwbcnt2")])
# visualizando estatísticas do modelo
summary(irt)

# gráfico de densidade dos valores dos itens
plot(factor.scores(irt), main = '', xlab = expression(theta), ylab = 'Densidade',
     lwd = 2)

# gráficos
plot(irt, type = 'ICC', lwd = 2, xlab = expression(theta), 
     ylab = 'Probabilidade', main = '', legend = T, cx = 'top', ncol = 5, cex.lab = 1.5)

# salvando fator no banco de dados
ThetaTable <- factor.scores(irt, 
                            resp.patterns = data.frame(lapply(ess[, c("imsmetn", "imdfetn", "impcntr", "imbgeco2", "imueclt2", 
                                                                      "imwbcnt2")], as.numeric)))[["score.dat"]]

ess <- cbind(ess, ThetaTable[, 8 : 10])
summary(ess$z1)

# scree plot
parallel <- ess %>% 
        dplyr::select(imsmetn, imdfetn, impcntr, imbgeco2, imueclt2, imwbcnt2) %>%
        fa.parallel(fm = 'minres', fa = 'fa', ylabel = 'Eigenvalue', main = '')

# fator e diagrama de correlação
factor <- ess %>% 
        dplyr::select(imsmetn, imdfetn, impcntr, imbgeco2, imueclt2, imwbcnt2) %>%
        fa(nfactors = 1, rotate = 'varimax')

# modificando o nome da variável latente
colnames(factor$loadings) <- 'Índice de 
Xenofobia'

# diagrama de correlação
fa.diagram(factor, main = '', digits = 2, e.size = .05, marg = c(.2, 5, .2, .2), 
           rsize = .4, labels = )

# criando banco de dados temporário para estimar correlações entre as variáveis 
# do índice
corr <- ess %>% dplyr::select(imsmetn, imdfetn, impcntr, imbgeco2, imueclt2, 
                              imwbcnt2) %>% 
        cor(method = 'spearman') %>% as.data.frame()

# gráfico de calor
heat_plot <- ggcorrplot(corr, outline.color = 'white', type = 'lower', hc.order = F)

# matriz de correlação
corr_mat <- ess %>% dplyr::select(imsmetn, imdfetn, impcntr, imbgeco2, imueclt2, 
                                  imwbcnt2) %>%
        ggpairs(lower = 'blank', 
                diag =  'blank', 
                upper = list(continuous = wrap("cor", method = "spearman"))) +
        theme_test()

ggsave('./Command files/plots/heat_plot.png', plot = heat_plot, device = 'png', width = 4, height = 3)
ggsave('./Command files/plots/corr_matrix.png', plot = corr_mat, device = 'png',
       width = 5, height = 4)

###### REGRESSÕES ######
# percepção sobre imigrantes
logit1 <- glmer(data = ess, vote_rrp ~ imwbcnt2 + (1 | cntry) + (1 | essround),
                family = 'binomial', weights = pspwght)
summary(logit1)

# número de imigrantes no país
logit2 <- glmer(data = ess, vote_rrp ~ noimbro + (1 | cntry) + (1 | essround),
                family = 'binomial', weights = pspwght)
summary(logit2)

# percepção sobre imigrantes + número de imigrantes no país
logit3 <- glmer(data = ess, vote_rrp ~ imwbcnt2 + noimbro + imwbcnt2*noimbro +
                        (1 | cntry) + (1 | essround), 
                family = 'binomial', weights = pspwght)
summary(logit3)

# percepção sobre imigrantes + número de imigrantes no país com pesos da 
# probabilidade de inserção do respondente na amostra
logit4 <- glmer(data = ess, vote_rrp ~ imwbcnt2 + noimbro + imwbcnt2*noimbro +
                        (1 | cntry) + (1 | essround), weights = dweight, family = 'binomial')
summary(logit4)

# índice de xenofobia
logit5 <- glmer(data = ess, vote_rrp ~ z1 + (1 | cntry) + (1 | essround), 
                family = 'binomial', weights = pspwght)
summary(logit5)

# interação entre o índice de xenofobia e número de imigrantes
logit6 <- glmer(data = ess, vote_rrp ~ z1 + noimbro + z1*noimbro + (1 | cntry) + 
                        (1 | essround), family = 'binomial', weights = pspwght)
summary(logit6)

# design weight
logit7 <- glmer(data = ess, vote_rrp ~ z1 + noimbro + z1*noimbro + (1 | cntry) + 
                        (1 | essround), family = 'binomial', weights = dweight)
summary(logit7)

# full models
# variável percepção sobre imigração
logit8 <- glmer(data = ess, vote_rrp ~ imwbcnt2 + noimbro + imwbcnt2*noimbro + 
                         stfdem + stfeco + stfgov + stflife + partisan + 
                         blue_c + gndr + agea + (1 | cntry) + (1 | essround), 
                 family = 'binomial', weights = pspwght)
summary(logit8)

# regressão com coeficientes padronizados
logit8_std <- glmer(data = ess, vote_rrp ~ scale(imwbcnt2) + scale(noimbro) + 
                            scale(imwbcnt2*noimbro) + scale(stfdem) + scale(stfeco) +
                            scale(stfgov) + scale(stflife) + scale(partisan) + 
                            scale(blue_c) + scale(gndr) + scale(agea) + (1 | cntry) + (1 | essround), 
                    family = 'binomial', weights = pspwght)
summary(logit8_std)

# índice de xenofobia
logit9 <- glmer(data = ess, vote_rrp ~ z1 + noimbro + z1*noimbro + stfdem + 
                        stfeco + stfgov + stflife + partisan + blue_c + gndr + 
                        agea + (1 | cntry) + (1 | essround), 
                family = 'binomial', weights = pspwght)
summary(logit9)

# logit 9 com coeficientes padronizados
logit9_std <- glmer(data = ess, vote_rrp ~ scale(z1) + scale(noimbro) + 
                            scale(z1*noimbro) + scale(stfdem) + scale(stfeco) +
                            scale(stfgov) + scale(stflife) + scale(partisan) + 
                            scale(blue_c) + scale(gndr) + scale(agea) + (1 | cntry) + (1 | essround), 
                    family = 'binomial', weights = pspwght)

summary(logit9_std)

# tabela com os coeficientes das regressões
stargazer(logit1, logit2, logit3, logit4, logit5, logit6, logit7, logit8, logit9, 
          out = 'coef.tex')

texreg(list(logit1, logit2, logit3, logit4, logit5, logit6, logit7, logit8, logit9),
       file = 'coef_texreg.tex', bold = 0.05)

###### GRÁFICOS DOS RESULTADOS ######
# transformando os modelos para formato "tidy"
l1_t <- tidy(logit1)
l2_t <- tidy(logit2)
l3_t <- tidy(logit3)
l4_t <- tidy(logit4)
l5_t <- tidy(logit5)
l6_t <- tidy(logit6)
l7_t <- tidy(logit7)
l8_t <- tidy(logit8)
l9_t <- tidy(logit9)
l8_std_t <- broom::tidy(logit8_std)
l9_std_t <- broom::tidy(logit9_std)


l1_t <- l1_t %>% mutate(model = '(1)') 
l2_t <- l2_t %>% mutate(model = '(2)') 
l3_t <- l3_t %>% mutate(model = '(3)') 
l4_t <- l4_t %>% mutate(model = '(4)') 
l5_t <- l5_t %>% mutate(model = '(5)') 
l6_t <- l6_t %>% mutate(model = '(6)') 
l7_t <- l7_t %>% mutate(model = '(7)') 
l8_t <- l8_t %>% mutate(model = '(8)') 
l9_t <- l9_t %>% mutate(model = '(9)')
l8_std_t <- l8_std_t %>% mutate(model = '(8)')
l9_std_t <- l9_std_t %>% mutate(model = '(9)')

# juntando modelos num único dataset
all_m <- bind_rows(l1_t, l2_t, l3_t, l4_t, l5_t, l6_t, l7_t, l8_t, l9_t)

# transformando coeficientes para probabilidades preditas
allm2 <- bind_rows(l1_t, l2_t, l3_t, l4_t, l5_t, l6_t, l7_t, l8_t, l9_t)
allm2$estimate <- exp(allm2$estimate) / (1 + exp(allm2$estimate))

# juntando modelos com coeficientes padronizados
allm_std <- bind_rows(l8_std_t, l9_std_t)

write.csv(allm2, './Analysis data/probmodels.csv')

# renomeando variáveis independentes
all_m <- all_m %>% relabel_predictors(c('z1' = 'Índice de Xenofobia',
                                        'z1:noimbro' = 'Índice de xenofobia*noimbro',
                                        'imwbcnt2' = 'imwbcnt', 
                                        'imwbcnt2:noimbro' = 'imwbcnt*noimbro'))
allm2 <- allm2 %>% relabel_predictors(c('z1' = 'Índice de Xenofobia',
                                        'z1:noimbro' = 'Índice de xenofobia*noimbro',
                                        'imwbcnt2' = 'imwbcnt', 
                                        'imwbcnt2:noimbro' = 'imwbcnt*noimbro'))

allm_std <- allm_std %>% relabel_predictors(c('scale(z1)' = 'Índice de Xenofobia',
                                              'scale(z1 * noimbro)' = 'Índice de xenofobia*noimbro',
                                              'scale(imwbcnt2)' = 'imwbcnt', 
                                              'scale(imwbcnt2 * noimbro)' = 'imwbcnt*noimbro'))

allm_filt <- all_m[all_m$term %in% c('Índice de Xenofobia', 'imwbcnt', 'noimbro'),]
allm2_filt <- allm2[allm2$term %in% c('Índice de Xenofobia', 'imwbcnt', 'noimbro'),]
allm_std_filt <- allm_std[allm_std$term %in% c('Índice de Xenofobia', 'imwbcnt'),]

# gráfico - coeficientes em log odds
coef1 <- dwplot(allm_filt, show_intercept = F, 
                vline = geom_vline(xintercept = 0, colour = 'grey 60', 
                                   linetype = 2)) +
        theme_sjplot2() +
        scale_color_discrete(name = 'Modelo') +
        labs(x = 'log odds(vote_rrp = 1)')

# gráfico - probabilidades preditas
coef2 <- dwplot(allm2_filt, show_intercept = F, 
                vline = geom_vline(xintercept = .5, colour = 'grey 60', 
                                   linetype = 2)) +
        theme_sjplot2() +
        scale_color_discrete(name = 'Modelo') +
        labs(x = 'P(vote_rrp = 1)')

# gráfico - coeficientes padronizados 
coef3 <- dwplot(allm_std_filt, show_intercept = F, 
                vline = geom_vline(xintercept = 0, colour = 'grey 60', 
                                   linetype = 2)) +
        theme_sjplot2() +
        scale_color_discrete(name = 'Modelo') +
        labs(x = expression(sigma^2))

# comparando coeficientes entre modelos
coef_com <- secret_weapon(all_m, var = 'Índice de Xenofobia') +
        theme_sjplot2() +
        labs(x = 'log odds(vote_rrp = 1)', y = 'Modelo') +
        scale_color_discrete('Variável', labels = 'Índice de
Xenofobia')

# gráficos de interação
# interação percepção sobre imigrantes - valores minimo e máximo
prob1 <- plot_model(logit8, type = 'int', mdrt.values = 'minmax') +
        labs(x = 'imwbcnt', y = 'P(vote_rrp)', title = '') +
        theme_sjplot2() +
        scale_y_continuous(limits = c(0, .85))

# percepção sobre imigrantes - média mais ou menos um desvio padrão
prob2 <- plot_model(logit8, type = 'int', mdrt.values = 'meansd') +
        labs(x = 'imwbcnt', y = 'P(vote_rrp)', title = '') +
        theme_sjplot2()+
        scale_y_continuous(limits = c(0, .85))

# índice de xenofobia - minimo e máximo
prob3 <- plot_model(logit9, type = 'int', mdrt.values = 'minmax') +
        labs(x = 'Índice de Xenofobia', y = 'P(vote_rrp)', title = '') +
        theme_sjplot2()+
        scale_y_continuous(limits = c(0, .85))

# índice de xenofobia - media e mais ou menos um desvio padrão
prob4 <- plot_model(logit9, type = 'int', mdrt.values = 'meansd') +
        labs(x = 'Índice de Xenofobia', y = 'P(vote_rrp)', title = '') +
        theme_sjplot2()+
        scale_y_continuous(limits = c(0, .85))

# efeito marginal do índice de xenofobia
prob5 <- plot_model(logit5, terms = 'z1', type = 'pred') +
        labs(x = 'Índice de Xenofobia', y = 'P(vote_rrp)', title = '') +
        theme_sjplot2()+
        scale_y_continuous(limits = c(0, .85))

# efeito marginal da percepção sobre imigrantes
prob6 <- plot_model(logit1, terms = 'imwbcnt2', type = 'pred') +
        labs(x = 'imwbcnt', y = 'P(vote_rrp)', title = '') +
        theme_sjplot2()+
        scale_y_continuous(limits = c(0, .85))

# criando função para plotar a função sigmoidal do Índice de Xenofobia
sigmoid = function(x) {
        1 / (1 + exp(-x))
}

plot(ess$z1, sigmoid(ess$z1))

# plotando função sigmoidal
sig_z1 <- ggplot(data = ess, aes(x = z1, y = sigmoid(z1))) +
        geom_smooth(se = T) +
        labs(y = 'P(vote_rrp = 1)', x = 'Índice de Xenofobia') +
        theme_sjplot2() +
        scale_y_continuous(labels = seq(0, 1, 0.25))

# excluindo NA's da variável dependente
ess2 <- ess %>% drop_na(vote_rrp)

# criando data frame para plotar os valores preditos por país e por ano
logit_plot <- data.frame(x = ess$z1, y = fitted(logit5), País = ess$cntry,
                         Ano = ess$essround)

# plotando gráfico
fit_logit5 <- ggplot(data = logit_plot, aes(x = x, y = y, colour = factor(Ano))) +
                geom_smooth() +
                facet_wrap(~País) +
                theme_sjplot2() +
                labs(x = 'Índice de Xenofobia', y = 'P(vote_rrp)') +
                scale_color_discrete(name = 'Ano', labels = seq(2002, 2016, 2))

# criando data frame para plotar os resultados do modelo 9 (full model)
logit_plot2 <- ess %>% drop_na(vote_rrp, z1, noimbro, stfdem, stfeco, stfgov, stflife, 
                               partisan, blue_c, gndr, agea)

# adicionando valores preditos ao banco de dados
logit_plot2$y <- fitted(logit9)

# filtrando as variáveis de interesse só por garantia
logit_plot2 <- logit_plot2 %>% dplyr::select(y, z1, essround, cntry)

table(logit_plot2$essround)

# plotando gráfico (modelo 9)
fit_logit9 <- ggplot(logit_plot2, aes(x = z1, y = y, colour = factor(essround))) +
        geom_smooth() +
        facet_wrap(~cntry) +
        labs(x = 'Índice de Xenofobia', y = 'P(vote_rrp)') +
        scale_color_discrete(name = 'Ano', labels = c(2002, 2014)) +
        theme_sjplot2()

# salvando gráficos
ggsave('coef_logodds.png', plot = coef1, device = 'png', width = 5, height = 3.5)
ggsave('coef_probs.png', plot = coef2, device = 'png', width = 5, height = 3.5)
ggsave('std_coef.png', plot = coef3, device = 'png', width = 5, height = 3.5)
ggsave('int_immig_minmax.png', plot = prob1, device = 'png', width = 5, height = 3)
ggsave('int_immig_sd.png', plot = prob2, device = 'png', width = 5, height = 3)
ggsave('int_IX_minmax.png', plot = prob3, device = 'png', width = 5, height = 3)
ggsave('int_IX_sd.png', plot = prob4, device = 'png', width = 5, height = 3)
ggsave('eff_IX.png', plot = prob5, device = 'png', width = 4.5, height = 3)
ggsave('eff_immig.png', plot = prob6, device = 'png', width = 4.5, height = 3)
ggsave('sigmoid_IX.png', plot = sig_z1, device = 'png', width = 4.5, height = 3)
ggsave('eff_IX_bycountry_5.png', plot = fit_logit5, device = 'png', width = 7, height = 4)
ggsave('eff_IX_bycountry_9.png', plot = fit_logit9, device = 'png', width = 7, height = 4)
ggsave('coef_com.png', plot = coef_com, device = 'png', width = 4.5, height = 3)



b01 <- all_m[all_m$term == '(Intercept)' & all_m$model == '(5)',]
b01 <- b01$estimate
z11 <- all_m[all_m$term == 'z1' & all_m$model == '(5)',]
z11 <- z11$estimate
z1_range <- seq(from = min(ess$z1), to = max(ess$z1), by = .001)

a_logit <- b01 + z11*z1_range

summary
a_probs <- exp(a_logit) / (1 + exp(a_logit))

plot(z1_range, a_probs)

###### TESTES DE ROBUSTEZ ######
# relogit
rlog <- zelig(data = ess, vote_rrp ~ z1 + noimbro + z1*noimbro + stfdem + 
                      stfeco + stfgov + stflife + partisan + blue_c + gndr + 
                      factor(cntry) + factor(essround),
              model = 'relogit')
summary(rlog)

# criando banco de dados com os coeficientes do relogit para as variáveis de interesse
# e juntando com as estimativas do modelo 9
rlog_df <- data.frame(term = c('z1', 'z1:noimbro'), 
                      estimate = c(0.898008, -0.005543),
                      std.error = c(0.050944, 0.001580),
                      statistic = c(17.627, -3.508),
                      p.value = c(0.0000000000000002, 0.000452),
                      group = 'fixed', 
                      model = 'relogit') %>% rbind(l9_t)

# transformando variáveis em caractere e renomeando variáveis
rlog_df$term <- as.character(rlog_df$term)
rlog_df$model <- as.character(rlog_df$model)
rlog_df$group <- as.character(rlog_df$group)
rlog_df$term[rlog_df$term == 'z1'] <- 'Índice de Xenofobia'
rlog_df$term[rlog_df$term == 'z1:noimbro'] <- 'Índice de Xenofobia*noimbro'

# plotando efeitos comparados com o modelo final = índice de xenofobia
rlog_p <- secret_weapon(rlog_df, var = 'Índice de Xenofobia', 
                        vline = geom_vline(xintercept = 0, colour = 'grey 60', 
                                                                linetype = 2)) +
        theme_sjplot2() +
        labs(x = 'log odds(vote_rrp = 1)', y = 'Modelo') 
        
# interação
rlog_p2 <- secret_weapon(rlog_df, var = 'z1:noimbro', 
                         vline = geom_vline(xintercept = 0, colour = 'grey 60', 
                                            linetype = 2))

# plotando o efeito das duas variáveis de interesse
relogplot <- small_multiple(rlog_df[rlog_df$term %in% c('Índice de Xenofobia', 'Índice de Xenofobia*noimbro'),]) +
        geom_hline(yintercept = 0, colour = "grey60", linetype = 2) +
        theme_sjplot2() +
        theme(legend.position = 'none') +
        labs(x = 'Modelo', y = 'log odds(vote_rrp = 1)')

ggsave('comp_relogit_model9.png', plot = relogplot, device = 'png', width = 5, height = 4.5)

ess2 <- ess2 %>% drop_na(vote_rrp,  z1 , noimbro, stfdem , 
                                 stfeco , stfgov , stflife , partisan , blue_c , gndr , 
                                 agea)

# construindo data frame para plotar o erro do modelo final
fit_se_p <- tibble(fit = fitted(logit9), 
                   se = resid(logit9),
                   Theta = ess2$z1,
                   Thetanoimbro = ess2$z1*ess2$noimbro,
                   vote_rrp = ess2$vote_rrp)

# plot de valores preditos x erro
fit_se <- ggplot() +
        geom_point(data = fit_se_p, aes(x = fit, y = se), alpha = .07) +
        geom_smooth(data = fit_se_p %>% filter(vote_rrp == 0), 
                    aes(x = fit, y = se), colour = 'red') +
        geom_smooth(data = fit_se_p %>% filter(vote_rrp == 1), 
                    aes(x = fit, y = se)) +
        geom_smooth(data = fit_se_p, aes(x = fit, y = se), method = 'auto', 
                    colour = 'black') +
        geom_hline(yintercept = 0, linetype = 2, colour = 'grey 60') +
        theme_sjplot2() +
        labs(x = 'Valores preditos', y = 'Resíduos')


# resíduos x índice de xenofobia
z1_se <- ggplot(data = fit_se_p) +
        geom_point(aes(x = Theta, y = se), alpha = .07) +
        geom_smooth(aes(x = Theta, y = se), formula = y ~ x, method = 'auto',
                    colour = 'black') +
        geom_smooth(data = fit_se_p %>% filter(vote_rrp == 1),
                    aes(x = Theta, y = se), formula = y ~ x) +
        geom_smooth(data = fit_se_p %>% filter(vote_rrp == 0),
                    aes(x = Theta, y = se), formula = y ~ x, colour = 'red') +
        geom_hline(yintercept = 0, linetype = 2, colour = 'grey 60') +
        theme_sjplot2() +
        labs(x = 'Índice de Xenofobia', y = 'Resíduos')

# resíduos vs. termo interativo
int_se <- ggplot(data = fit_se_p) +
        geom_point(aes(x = Thetanoimbro, y = se), alpha = .07) +
        geom_smooth(aes(x = Thetanoimbro, y = se), formula = y ~ x, method = 'auto')+
        geom_hline(yintercept = 0, linetype = 2, colour = 'grey 60') +
        geom_vline(xintercept = 0, linetype = 2, colour = 'grey 60') +
        theme_sjplot2() +
        labs(x = 'Índice de Xenofobia*noimbro', y = 'Resíduos')

ggsave('fitse.png', plot = fit_se, device = 'png', width = 5, height = 3.5)
ggsave('z1se.png', plot = z1_se, device = 'png', width = 5, height = 3.5)
ggsave('z1_noimbrose.png', plot = int_se, device = 'png', width = 5, height = 3.5)

# binnedplot modelo 9
binnedplot(fitted(logit9), resid(logit9, type = 'response'),
           xlab = 'Valores preditos', ylab = 'Resíduo médio', main = '',
           cex.pts = .8, col.pts = 1, col.int = 'gray')

# computando vif
vif_9 <- vif(logit9) %>% as.tibble()
vif_9$Variável <- c('Índice de Xenofobia', 'noimbro', 'stfdem', 'stfeco', 'stfgov',
                    'stflife', 'partisan', 'blue_c', 'gender', 'agea', 'Interação')

# plotando num gráfico de barra
bar_vif <- ggplot(data = vif_9, aes(x = reorder(Variável, -value), y = value)) +
        geom_bar(stat = 'identity') +
        theme_sjplot2() +
        geom_hline(yintercept = 5, linetype = 2, colour = 'grey 60') +
        labs(y = 'VIF', x = '') +
        theme(axis.text.x = element_text(angle = 90))

ggsave('vif.png', plot = bar_vif, device = 'png', width = 4.5, height = 3)
?vif

# computando residuos padronizados e cook's distance
model.data <- augment(logit9) %>% mutate(index = 1:n())

# function to name the top 3 largest values
model.data %>% top_n(3, .cooksd)

# ploting std. residuals
stdresid <- ggplot(model.data, aes(index, .std.resid)) + 
        geom_point(aes(color = vote_fn), alpha = .5) +
        theme_classic() +
        ggtitle("Observations vs. Standardized residuals - Model 4") +
        labs(x = "Index", y = "Standardized residuals", 
             caption = "Figure 8 - Author's elaboration") + 
        theme(plot.title = element_text(size = 12, hjust = 0.5, face = "bold"),
              legend.title = element_text(size = 10, face = "bold"),
              plot.caption = element_text(hjust = 0.5),
              legend.position = "right",
              legend.title.align = .5) + 
        scale_color_discrete(name = "Vote in the FN")